import xarray as xr
import wrf
import netCDF4 as nc
import os 
import numpy as np
import pandas as pd
from scipy.interpolate import griddata
import metpy
from metpy.units import units 
import sys



import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker  # 导入刻度格式化器
import matplotlib.colors as mcolors


### settings ###
wrfpath='/public23/home/sc61341/gsj/wrfout/wrfout_d01_2023-09-06_00:00:00'
fmpath='/public23/home/sc61341/gsj/FM_1.0/test/em_real/wrfout_d01_2023-07-29_00:00:00'
obsdirname='/public23/home/sc61341/gsj/data/23BJmetstation/r24_MetStation'  #不需要给最后的/
era5path='/public23/home/sc61341/gsj/data/ERA5/case/beijing729/beijing729.nc'
plotlevels=[500,700,850]  #需要绘制的垂直高度层次；
casename="beijing_newter"  # 创建casename文件放图
cat='/'  #windows--\\  linxu--/ 

########################################

path=os.path.dirname(__file__)
def makefiletree(path,casename):
    path+=f"{cat}"+casename
    if not os.path.exists(path):
        os.mkdir(path)
    model=[f'{cat}WRF',f'{cat}FM',f'{cat}ERA5',f'{cat}OBS']
    fig=[f'{cat}H_T_UV',f'{cat}W_vaportrans',f'{cat}rainsnow']
    for m in model:
        dirname=path+m
        if not os.path.exists(dirname):
            os.mkdir(dirname)       
        for f in fig:
            dirname=path+m+f
            if not os.path.exists(dirname):
                os.mkdir(dirname)        
makefiletree(path,casename)

class metplot:
    def __init__(self):
        fig, ax = plt.subplots(subplot_kw={'projection': ccrs.PlateCarree()},figsize=(15, 10),dpi=400)
        self.fig=fig
        self.ax=ax
        #rainsetting
        self.rain_levels=[0.1,10.0,25.0,50.0,100.0,250.0,500.0]#雨量等级
        raindict=['#A6F28F','#3DBA3D','#61BBFF','#0000FF','#FA00FA','#800040']#颜色列表
        self.rain_cmap=mcolors.ListedColormap(raindict)#产生颜色映射
        self.rain_norm=mcolors.BoundaryNorm(self.rain_levels,self.rain_cmap.N)#生成索引
        #snow setting
        self.snow_levels = [0, 0.1, 2.5, 5, 10, 25, 50] #雪量
        snowcolors = ["#ffffff", "#c5e7b0", "#68b165", "#329ed8", "#286ab3", "#ec46f6"]
        self.snow_cmap = mcolors.ListedColormap(snowcolors)
        self.snow_norm=mcolors.BoundaryNorm(self.snow_levels,self.snow_cmap.N)#生成索引
        
        #底图；lon,lat,title代表底图范围和标题（画图函数可以比底图大，但只显示地图）
    def base(self,lon,lat,title=None):
        #title+map
        self.ax.add_feature(cfeature.COASTLINE) #海岸线
        self.ax.add_feature(cfeature.BORDERS, linestyle=':')   #国界
        china_provinces = cfeature.NaturalEarthFeature(category='cultural', name='admin_1_states_provinces_lines', scale='10m', facecolor='none')
        self.ax.add_feature(china_provinces, edgecolor='black', linewidth=0.5) #省界
        self.ax.set(xlim=(np.min(lon),np.max(lon)))
        self.ax.set(ylim=(np.min(lat),np.max(lat)))    
        self.ax.set_xticks(np.linspace(np.min(lon), np.max(lon), 6))  
        self.ax.xaxis.set_major_formatter(cticker.LongitudeFormatter(number_format='.1f'))  # 格式化为经度格式
        self.ax.set_yticks(np.linspace(np.min(lat), np.max(lat), 4))  
        self.ax.yaxis.set_major_formatter(cticker.LatitudeFormatter(number_format='.1f'))  # 格式化为纬度格式
        self.ax.set_title(title)

        
    def contourf(self,lon,lat,ctrf,cmap="bwr",levels=None,norm=None):
        contourf=self.ax.contourf(lon,lat,ctrf,cmap=cmap,alpha=1,levels=levels,extend="both",norm=norm)
        self.fig.colorbar(contourf, ax=self.ax, shrink=0.7, orientation='vertical' )
    
    def contour(self,lon,lat,ctr,levels=None,norm=None,colors='black'):
        contour = self.ax.contour(lon, lat , ctr, levels=levels, norm=norm, colors=colors)  # 绘制等高线图
        self.ax.clabel(contour, inline=True, fontsize=10)  # 为等高线添加标签
    
    def barb(self,lon,lat,u,v,barbcolor=None):
        sizes = dict(emptybarb=0.1, spacing=0.15, height=0.6)  #height 外凸长度；empytbarb0风速大小；spacing间隔；值均是相对于length的比例
        self.ax.barbs(lon, lat, u, v, transform=ccrs.PlateCarree(),
              sizes=sizes, length=6,barb_increments={'half':2,'full':4,'flag':20}, #length长度
              linewidth=1.6, regrid_shape=10,barbcolor=barbcolor)
    # WRF点太密，对barb quiver图做下采样    
    def sample_factor(self,lon,lat):
        n=lon.ndim
        match n:
            case 2:
                ny,nx=lon.shape
                ns=np.min([ny,nx])/10
            case 1:
                nx=lon.shape[0]
                ny=lat.shape[0]
                ns=np.min([ny,nx])/10
            case _:
                print('!!!!!! new situation appears!!!!! need to update sampling function !!!!!')
                os._exit(0)
        return int(ns)
        
    def quiver(self,lon,lat,u,v,color='green',ref=0.01,unit=None,scale=5e-5,width=0.003):
        n=self.sample_factor(lon, lat)
        match lon.ndim:
            case 2:        
                quiver=self.ax.quiver(lon[::n,::n],lat[::n,::n], u[::n,::n], v[::n,::n],color=color,scale=scale,scale_units='dots',width=width)
            case 1:
                n=3
                quiver=self.ax.quiver(lon[::n],lat[::n], u[::n,::n], v[::n,::n],color=color,scale=scale,scale_units='dots',width=width)
        self.ax.quiverkey(quiver,                      #传入quiver句柄
          X=0.9, Y = 0.9,       #确定 label 所在位置，都限制在[0,1]之间
          U = ref,                    #参考箭头长度 
          angle = 0,            #参考箭头摆放角度。默认为0，即水平摆放
          label=f'{ref} {unit}',        #箭头的补充：label的内容  + 
          labelpos='S',          #label在参考箭头的哪个方向; S表示南边
          color = color,labelcolor = 'b', #箭头颜色 + label的颜色
          ) 

    def save(self,path,title):
        plt.savefig(f'{path}{cat}{title}',dpi=400, bbox_inches='tight')
        plt.close()
        
    def glance(self,var):
        os.system("mkdir glance") ; os.chdir('glance')
        if var.ndim==2:
            plt.contourf(var);plt.colorbar()
            plt.savefig("1")
            plt.clf()
        elif var.ndim==3:
            for p in range(var.shape[0]):
                plt.contourf(var[p,:,:]);plt.colorbar()
                plt.savefig(f"{p}")
                plt.clf()    
        
        

class WRFCLASS(metplot):
    def __init__(self,wrfpath,plotlevels):
        self.f=nc.Dataset(wrfpath)
        self.Times=[f[:13] for f in wrf.getvar(self.f,"Times",timeidx=wrf.ALL_TIMES).values.astype(str)]
        self.plotlevels=plotlevels
        self.lon=wrf.getvar(self.f,"XLONG")
        self.lat=wrf.getvar(self.f,"XLAT")
        
        
    def overview(self):
        print(self.f.variables.keys())
        print("**********************")
        print("**********************")
        print(self.f["PH"])
        print("**********************")
        print("**********************")
        print(" destagger_usage :  wrf.destagger( wrf.getvar(f,'W').values,0) ")
        print(" getvar_usage :  wrf.destagger( p=wrf.getvar(f,'pressure',timeidx=t)  ps: timeidx=wrf.ALL_TIMES")
        print("interplevel_usage : wrf.interplevel(V,p,level)  ")
        print("**********************")
        print("**********************")
    
    def itpsmth(self,strvarname,itplevel,t):
        from scipy.ndimage import uniform_filter
        P=wrf.getvar(self.f,"pressure",timeidx=t)
        var=wrf.getvar(self.f,strvarname,timeidx=t)
        var_itp=wrf.interplevel(var,P,itplevel)
        match bool(np.isnan(var_itp).values.any() ):
            case True:
                var_itp_smooth=var_itp
            case False:
                var_itp_smooth=uniform_filter( var_itp, size=6)
        return var_itp_smooth
  
    def fig1(self,trange):
        # 1)hgt+temperature(ctr-black-red); uv(jet特标); ---温度平流
        lon=wrf.getvar(self.f,"XLONG")
        lat=wrf.getvar(self.f,"XLAT")
        for t in trange:
            for p in self.plotlevels:
                super().__init__()
                title=f"WRF___{p}hpa_{self.Times[t]}_Temp(redctr)_geopt(bluectr)_uv(barb)"
                geopt=self.itpsmth("geopt", p, t)
                T=self.itpsmth("tk", p, t)-273.15
                u=self.itpsmth("ua",p,t)
                v=self.itpsmth("va",p,t)
                self.base(lon,lat,title)
                self.contour(lon,lat,geopt,colors='blue')
                self.contour(lon,lat,T,colors="red")
                self.barb(lon,lat,u,v)
                self.save(path+cat+casename+cat+"WRF"+cat+'H_T_UV',title)
                
    def fig2(self,tarange):
        # 2) w+s水汽输送
        lon=wrf.getvar(self.f,"XLONG")
        lat=wrf.getvar(self.f,"XLAT")
        for t in tarange:
            for p in self.plotlevels:
                super().__init__()
                title=f"WRF___{p}hpa_{self.Times[t]}_W(contourf)_quiver(vaport_transport))"
                W=self.itpsmth("wa", p, t)           
                w=self.itpsmth("QVAPOR",p,t)
                q=w/(1+w) #混合比 to 比湿
                Uw=self.itpsmth("ua",p,t)*q/9.8
                Vw=self.itpsmth("va",p,t)*q/9.8
                                     
                self.base(lon,lat,title)
                self.contourf(lon,lat,W,levels=np.linspace(-2.4,2.4,9),cmap="seismic")
                self.quiver(lon,lat,Uw,Vw,ref=0.01,unit="kg/m•hPa•s",scale=5e-5)
                self.save(path+cat+casename+cat+"WRF"+cat+'W_vaportrans',title)
            
    def fig3(self,trange):
        # 24h雨/雪 trange 只给出开始时刻
        lon=wrf.getvar(self.f,"XLONG")
        lat=wrf.getvar(self.f,"XLAT")
        rain=( wrf.getvar(self.f,"RAINC",timeidx=wrf.ALL_TIMES) + 
            wrf.getvar(self.f,"RAINNC",timeidx=wrf.ALL_TIMES) +
            wrf.getvar(self.f,"RAINSH",timeidx=wrf.ALL_TIMES) )
        snow=wrf.getvar(self.f,"SNOW",timeidx=wrf.ALL_TIMES)
        #rain
        match len(trange):
            case 1:
                print('Wrong when drawing snow/rain')
                print('Input trange is wrong,you should at least input a list containing start and strat+gap')
                sys.exit()
            case _:
                gap=trange[1]-trange[0]
        for t in trange:
            if t+gap>len(self.Times):
                pass
            else:
                super().__init__()
                title=f"WRF_Rain24h___{self.Times[t]}-{self.Times[t+gap]}"
                self.base(lon,lat,title)
                self.contourf(lon,lat,rain[t+gap]-rain[t],cmap=self.rain_cmap,norm=self.rain_norm,levels=self.rain_levels )
                self.save(path+cat+casename+cat+'WRF'+cat+'rainsnow',title)
        #snow
        for t in trange:
            if t+gap>len(self.Times):
                pass
            else:
                super().__init__()
                title=f"WRF_Snow24h___{self.Times[t]}-{self.Times[t+gap]}"
                self.base(lon,lat,title)
                self.contourf(lon,lat,snow[t+gap]-snow[t],cmap=self.snow_cmap,norm=self.snow_norm,levels=self.snow_levels )
                self.save(path+cat+casename+cat+'WRF'+cat+'rainsnow',title)


class OBSCLASS(metplot):
    #默认是24h一次的 txt文档； 只能画降水。
    def __init__(self,obsdirname,lon,lat): #lon,lat 要插值到的网格上(一般是WRF网格);
        self.fnames=os.listdir(obsdirname)
        self.fpaths=[obsdirname+cat+f for f in self.fnames]
        self.sep=0.2
        self.lon=lon
        self.lat=lat
    def interp24(self,filepath):
        t=pd.read_csv( filepath,encoding='gbk',skiprows=range(11),delim_whitespace=True,header=None)  #############
        lon_sta=t.iloc[:,1]   ############   
        lat_sta=t.iloc[:,2]
        rain_sta=t.iloc[:,4]
        A=np.column_stack((lon_sta,lat_sta))   ############  tuple form
        rain=griddata(A,rain_sta,(self.lon,self.lat),method="cubic")
        return rain
        
    def interp1h(self): #需要再补充
        pass 
    
    def rain(self):
        # 默认是所有天数的降水全画
        for t in range(len(self.fpaths)):
            super().__init__()
            filepath=self.fpaths[t]
            rain=self.interp24(filepath)
            title=f"OBS-24hRAIN_before_{self.fnames[t][:8]}"
            self.base(self.lon,self.lat,title)
            self.contourf(self.lon,self.lat,rain,cmap=self.rain_cmap,norm=self.rain_norm,levels=self.rain_levels )
            self.save(path+cat+casename+cat+'OBS'+cat+'rainsnow',title)
        
    def snow(self):
        for t in range(len(self.fpaths)): 
            super().__init__()
            filepath=self.fpaths(t)
            snow=self.interp24(filepath)
            title=f"OBS-24hSNOW_before_{self.fnames[t][:8]}"
            self.base(self.lon,self.lat,title)
            self.contourf(self.lon,self.lat,snow,cmap=self.snow_cmap,norm=self.snow_norm,levels=self.snow_levels )
            self.save(path+cat+casename+cat+'OBS'+cat+'rainsnow',title)
            

class ERA5CLASS(metplot):
    
    def __init__(self,era5path,lon,lat,plotlevels):
        self.ds=xr.open_dataset(era5path, engine="netcdf4")
        self.wrflon=lon
        self.wrflat=lat
        self.lon=self.ds.longitude.values
        self.lat=self.ds.latitude.values
        self.P=self.ds.pressure_level.values
        self.time=[str(t)[:13]  for t in self.ds.valid_time.values]
        self.plotlevels=plotlevels
    
    def overview(self):
        print(self.ds)
        print(self.time)
        
    def get2d(self,strvarname,t,p):
        return self.ds[strvarname][t,:,:,:].loc[p,:,:]
    
    
    def fig1(self,trange):
        for t in trange:
            for p in self.plotlevels:
                super().__init__()
                T=self.get2d("t", t, p)-273.15
                H=self.get2d("z",t,p)
                U=self.get2d('u',t,p)
                V=self.get2d('v',t,p)             
                title=f"{self.time[t]}_{p}hPa_T(contour-red)_H(contour-black)_UV(barb) "
                self.base(self.wrflon,self.wrflat,title)
                self.contour(self.lon,self.lat,H,colors='blue')
                self.contour(self.lon,self.lat,T,colors="red")
                self.barb(self.lon,self.lat,U,V)
                self.save(path+cat+casename+cat+"ERA5"+cat+'H_T_UV',title)
    
    def fig2(self,trange):
        # 静力平衡假设
        for t in trange:
            for p in self.plotlevels:
                super().__init__()
                title=f"ERA5_{self.time[t]}_{p}hPa_W(contourf)_quiver(vaport_transport)"
                # calc W
                omega=self.get2d("w",t,p)*units('Pa/s')
                T=self.get2d("t", t, p)*units('K')
                q=self.get2d('q', t, p)*units('kg/kg')
                W=metpy.calc.vertical_velocity(omega,p*units('hPa'),T,q)
                #calc 
                Uq=self.get2d('u',t,p)*q/9.8
                Vq=self.get2d('v',t,p)*q/9.8
                self.base(self.wrflon,self.wrflat,title)
                self.contourf(self.lon,self.lat,W,levels=np.linspace(-2.4,2.4,9),cmap="seismic")
                self.quiver(self.lon, self.lat, Uq, Vq,ref=0.01,unit="kg/m•hPa•s",scale=5e-5)
                self.save(path+cat+casename+cat+"ERA5"+cat+"W_vaportrans",title)
        
class FMCLASS(metplot):
    def __init__(self,fmpath,plotlevels):
        self.f=nc.Dataset(fmpath)
        self.Times=[f[:13] for f in wrf.getvar(self.f,"Times",timeidx=wrf.ALL_TIMES).values.astype(str)]
        self.lon=wrf.getvar(self.f,"XLONG").values
        self.lat=wrf.getvar(self.f,'XLAT').values
        self.z=wrf.destagger( wrf.getvar(self.f,'ZFULL'),0)
        self.plotlevels=plotlevels
        
    def itp(self,var,t,p):
        P=wrf.getvar(self.f,'PFM',timeidx=t)/100
        return wrf.interplevel(var,P,p)    
    
    def smooth2d(self,var):
        from scipy.ndimage import uniform_filter
        match bool(np.isnan(var).values.any()): #np.isnan返回numpy的bool_类型;这个类型不会被内置的python match case语句识别
            case True:
                return var
            case False:
                var=uniform_filter( var, size=6)
                return var

    def fig1(self,trange):
        # 1)hgt+temperature(ctr-black-red); uv(jet特标); ---温度平流
        for t in trange:
            for p in self.plotlevels:
                super().__init__()
                title=f"FM___{p}hpa_{self.Times[t]}_Temp(redctr)_geopt(bluectr)_uv(barb)"
                #get all var
                P=wrf.getvar(self.f,'PFM',timeidx=t)
                RHO=wrf.getvar(self.f,'RHOFM',timeidx=t)
                w=wrf.getvar(self.f,"QVAPOR",timeidx=t)
                q=w/(1+w)
                T=P/RHO/287.04/(1+0.618*q)
                U=wrf.getvar(self.f,'ua',timeidx=t)
                V=wrf.getvar(self.f,'va',timeidx=t)
                z=self.z
                # interpolation
                Tp=self.itp(T,t,p)
                Up=self.itp(U,t,p)
                Vp=self.itp(V,t,p)
                geopt=self.itp(z,t,p)/9.8
                # smooth
                Tp=self.smooth2d(Tp)
                geopt=self.smooth2d(geopt)
                # plot
                self.base(self.lon,self.lat,title)
                self.contour(self.lon,self.lat,geopt,colors='blue')
                self.contour(self.lon,self.lat,Tp,colors="red")
                self.barb(self.lon,self.lat,Up,Vp)
                self.save(path+cat+casename+cat+"FM"+cat+'H_T_UV',title)
                
    def fig2(self,trange):
        # 2) w+s水汽输送
        for t in trange:
            for p in self.plotlevels:
                super().__init__()
                title=f"FM___{p}hpa_{self.Times[t]}_W(contourf)_quiver(vaport_transport))"
                #getvar
                W=wrf.getvar(self.f,"wa", timeidx=t)           
                w=wrf.getvar(self.f,"QVAPOR",timeidx=t)
                U=wrf.getvar(self.f,'ua',timeidx=t)
                V=wrf.getvar(self.f,'va',timeidx=t)
                #interplevel
                W=self.itp(W,t,p) #这里的p和T的计算有关。不是因为W已取特定时刻而是多余的。
                w=self.itp(w,t,p)
                q=w/(1+w) #混合比 to 比湿
                U=self.itp(U,t,p)
                V=self.itp(V,t,p)
                # vapor transport
                Uw=U*q/9.8
                Vw=V*q/9.8
                # draw
                self.base(self.lon,self.lat,title)
                self.contourf(self.lon,self.lat,W,levels=np.linspace(-2.4,2.4,9),cmap="seismic")
                self.quiver(self.lon,self.lat,Uw,Vw,ref=0.01,unit="kg/m•hPa•s",scale=5e-5)
                self.save(path+cat+casename+cat+"FM"+cat+'W_vaportrans',title)
            
    def fig3(self,trange):
        # 24h雨/雪 trange 只给出开始时刻
        rain=( wrf.getvar(self.f,"RAINC",timeidx=wrf.ALL_TIMES) + 
            wrf.getvar(self.f,"RAINNC",timeidx=wrf.ALL_TIMES) +
            wrf.getvar(self.f,"RAINSH",timeidx=wrf.ALL_TIMES) )
        snow=wrf.getvar(self.f,"SNOW",timeidx=wrf.ALL_TIMES)
        #rain
        match len(trange):
            case 1:
                print('Wrong when drawing snow/rain')
                print('Input trange is wrong,you should at least input a list containing start and strat+gap')
                sys.exit()
            case _:
                gap=trange[1]-trange[0]
        for t in trange:
            if t+gap>len(self.Times):
                pass
            else:
                super().__init__()
                title=f"FM_Rain24h___{self.Times[t]}-{self.Times[t+gap]}"
                self.base(self.lon,self.lat,title)
                self.contourf(self.lon,self.lat,rain[t+gap]-rain[t],cmap=self.rain_cmap,norm=self.rain_norm,levels=self.rain_levels )
                self.save(path+cat+casename+cat+'FM'+cat+'rainsnow',title)
        #snow
        for t in trange:
            if t+gap>len(self.Times):
                pass
            else:
                super().__init__()
                title=f"FM_Snow24h___{self.Times[t]}-{self.Times[t+gap]}"
                self.base(self.lon,self.lat,title)
                self.contourf(self.lon,self.lat,snow[t+gap]-snow[t],cmap=self.snow_cmap,norm=self.snow_norm,levels=self.snow_levels )
                self.save(path+cat+casename+cat+'FM'+cat+'rainsnow',title)
                        


'''
#### Usage guide #########

Short guide : 看源码

Long guide:    
    为可视化模式性能 ；绘制同一个例的 ERA5/FM/WRF/OBS_station的天气形势和风场（FIG1）；垂直速度和水汽通量（FIG2）；降水/雪量（FIG3）。
    使用前，首先在代码开头设置一些参数；四种数据目录，绘制高空图的气压层；设置linux或windows环境的路径连接符
    画图代码简介：
        OBS:        
            OBS.rain()或OBS.snow() 不用传时间参数，直接绘制所有的。但是具体代码需要根据手头的txt资料改动。
        WRF/FM:
            自己传入需要绘制的时刻; 时刻通常以np.arange形式给出；
            先使用WRF.Times和FM.Times看时次。然后自己算一下自己想要哪个时次；
                其中WRF.fig3和FM.fig3中的时次是指的开始时间；会绘制开始时间后24h内的降水
            使用例子：
                绘制单个时刻测试画图效果  WRF.fig1( [1] )
                每隔6h绘制一次  WRF.fig1(np.arange(start,end,6) )
                每个24h绘制一次降水/雪 WRF.fig3(np.arange(start,end,24))
        ERA5：
            需要传绘制时刻; 但是由于ERA5一般是下载6h一次或者3h一次的。可以直接全画。 ex： ERA5.fig1(range(len(ERA5.time)))
    
    全时刻绘制代码： 若WRF,FM-1h ；ERA-6h ； OBS-24h
        WRF.fig1(np.arange(0,len(WRF.Times),6))
        WRF.fig2(np.arange(0,len(WRF.Times),6))
        WRF.fig3(np.arange(0,len(WRF.Times),24))
        FM.fig1(np.arange(0,len(FM.Times),6))
        FM.fig2(np.arange(0,len(FM.Times),6))
        FM.fig3(np.arange(0,len(FM.Times),24))
        ERA5.fig1(range(len(ERA5.time)))
        ERA5.fig2(range(len(ERA5.time)))
        OBS.rain()
    
'''
WRF=WRFCLASS(wrfpath,plotlevels)
ERA5=ERA5CLASS(era5path,WRF.lon,WRF.lat,plotlevels)
OBS=OBSCLASS(obsdirname,WRF.lon,WRF.lat)
FM=FMCLASS(fmpath,plotlevels)

###DIY wanted pic  #######
'''
WRF.fig1(np.arange(0,len(WRF.Times),1))
WRF.fig2(np.arange(0,len(WRF.Times),1))
WRF.fig3(np.arange(0,len(WRF.Times),4))
FM.fig1(np.arange(0,len(FM.Times),1))
FM.fig2(np.arange(0,len(FM.Times),1))
FM.fig3(np.arange(0,len(FM.Times),4))
ERA5.fig1(range(len(ERA5.time)))
ERA5.fig2(range(len(ERA5.time)))
OBS.rain()
'''

FM.fig1(np.arange(0,len(FM.Times),1))
FM.fig2(np.arange(0,len(FM.Times),1))
FM.fig3(np.arange(0,len(FM.Times),4))
